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ABSTRACT 

Both ground and space-based transit observatories are poised to significantly 
increase the number of known transiting planets and the number of precisely 
measured transit times. The variation in a planet's transit times may be used 
to infer the presence of additional planets. Deducing the masses and orbital 
parameters of such planets from transit time variations (TTVs) alone is a rich 
and increasingly relevant dynamical problem. In this work, we evaluate the ex- 
tent of the degeneracies in this process, systematically explore the dependence 
of TTV signals on several parameters and provide phase space plots that could 
aid observers in planning future observations. Our explorations are focused on 
a likely-to-be prevalent situation: a known transiting short-period Neptune or 
Jupiter-sized planet and a suspected external low-mass perturber on a nearly- 
coplanar orbit. Through ~ 10^ N-body simulations, we demonstrate how TTV 
signal amplitudes may vary by orders of magnitude due to slight variations in any 
one orbital parameter (10~^ AU in semimajor axis, 0.005 in eccentricity, or a few 
degrees in orbital angles) , and quantify the number of consecutive transit obser- 
vations necessary in order to obtain a reasonable opportunity to characterize the 
unseen planet (> 50 observations). Planets in or near period commensurabilities 
of the form p:q, where p < 20 and g < 3, produce distinct TTV signatures, re- 
gardless of whether the planets are actually locked in a mean motion resonance. 
We distinguish these systems from the secular systems in our explorations. Ad- 
ditionally, we find that computing the autocorrelation function of a TTV signal 
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can provide a useful diagnostic for identifying possible orbits for additional plan- 
ets and suggest that this method could aid integration of TTV signals in future 
studies of particular exosystems. 

Subject headings: celestial mechanics — methods: n-body simulations, statistical 
— stars: planetary systems 



Introduction 



1.1. Motivation 



Of the first ~ 60 exoplanets detected by transit photometry, none were accompanied by 



additional planets in the same system. The discovery of HAT-P-13c (Bakos et al. 2009) and 



CoRoT-7 dQueloz et al?][2009| broke new ground as the first systems to contain both a planet 
observed to transit and a second planet detected by other means. Then Kepler-9 became the 
first (and so far only) confirmed system with multiple transiting planets ( Holman et al.|20T0 ), 
although five other Kepler target stars now show evidence of multiple transiting candidates 
( Staff en et al.|[2010| ). With nearly a third of all known transiting planets formally announced 
in the year 200S|3 the relentless pace of transit detections suggests that soon investigators 
will discover an abundance of multi-planet systems containing at least one known transiting 
planet. 

This trend in exoplanet astrophysics highlights the importance of performing follow-up 
observations for 1) single-planet transiting systems, as a way to determine if additional plan- 
ets exist, and 2) multi-planet systems with at least one transiting planet, in order to better 
constrain the parameters of all planets in those systems. The gravitational tug of planetary 
perturbers on a known transiting Hot Jupiter will cause variations in the mid-transit times of 
that planet. Encoded in these transit timing variations (TTVs) is the influence of the hidden 
planet (s), and the process of extracting the mass and orbital parameters of these planets will 
become increasingly important as the number of transit detections increases. Up to now, the 
only planets definitiely confirmed by TTVs is for a system where both planets transit their 



parent star (Holman et al. 2010). 
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1.2. Observational Studies 



Previous studies have attempted to solve tliis "inverse" problem of deducing the mass 
and orbital parameters of an unseen perturber from a limited number of TTVs for a few spe- 



cific systems. Steffen & Agol (2005) performed one of the first TTV analysis on a particular 
exosystem when they studied 12 transit observations for TrES-1, and demonstrated that the 
data could have identified a hypothetical perturber in that system that is at the order of an 



Earth-mass, or lower. Agol & Steffen (2007) then combined 13 transit observations with 68 



radial velocity measurements for HD 209458 in order to constrain the presence of additional 



planets in that system. Miller- Ricci et al. (2008a) obtained up to 12 consecutive transits 



from MOST data sets of HD 209458b in 2004 and 2005, and obtained no TTV signatures 



above 80 s. Miller- Ricci et al. (2008b) obtained up to 10 consecutive transits from MOST 
data sets of HD 189733b in 2006, and obtained no TTVs above 45 s. The authors of both 
papers then used these parameters in order to restrict the possible existence of additional 



planets in these systems. Diaz et al. (2008) analyzed just 5 transit data points to conclude 
that OGLE-TR-lllb cannot produce the variations in signal; they rule out the presence of 
a satellite causing the TTV, and instead suggest an exterior Earth-mass planet could be 



the source of the variation. Adams et al. (2010), using data from 6 additional tranits, has 



subsequently challenged these claims. Coughlin et al. (2008) analyzed 28 transit observa- 



tions for the approximately Neptune-mass planet, Gl 436b, and were able to rule out the 



existence of any planets which cause a TTV of over 60 s. Csizmadia et al. (2010) considered 
the TTV on 36 transits of CoRoT-lb and did not find any periodic signals, so as to rule out 
additional planets in the form of Super Earths, Saturn-like planets and Jupiter-like planets 
each in particular regions of parameter space. 

All these investigations report fewer than 50 transits. We will demonstrate that in 
most cases, at least 50 transits are needed in order to appreciably narrow the phase space 
of possible solutions to the inverse problem. However, the ongoing space-based missions 
CoRoT and Kepler provide cause for optimism. CoRoT has already discovered at least 7 



planets ( 


Barge et al.||2008 


Alonso et al. 


2008 




Aigrain et al. 


2008 


Deleuil et al. 


2008 


Rauer 


et al. 


2009 


Leger et al. 


2009 


Fridlund et al. 


2010 


) and Kepler at least 5 ( 


Borucki et al. 



2010). The nearly continuous observational coverage for 0.5 — 3.5 yr by both missions should 



allow for hundreds of transits of a single planet to be observed. 



1.3. Theoretical Studies 



Purely theoretical studies have appealed to both analytical considerations and N-body 
simulations. The nature of the inverse problem coupled with the observed constraints on the 



-4- 



architecture of observable planetary systems with a transiting planet dictate which relevant 
branches of perturbation theory are useful for TTV studies. When N-body simulations are 
used, they need to be run for half of one year to several years. Such constraints help provide 
the context for the several important theoretical contributions which have established a 
foundation on which future TTV studies may be based and which we now proceed to outline. 



Holman & Murray (2005) and Agol et al. (2005), in two early studies, considered sys 



tems with two coplanar planets. Holman & Murray (2005) reported a) that disorderly and 
aperiodic TTV signals may arise from the presence of multiple planets (their Fig. 1), b) on 
the correlation between planetary period ratio and timing signal for fixed values of several 
other orbital parameters and masses (their Figs. 3-4), and c) an estimate for the amplitude 
of this signal using Laplace-Lagrange secular theory (their Eq. 1). These results helped 
demonstrate that secular theory can produce the correct order of magnitude signal in many 
cases, but fail to describe the spikes in TTV signal which appear when the planetary period 
ratio is commensurate. 



Agol et al. (2005) explored a greater region of phase space, and focused their investiga- 



tions on four regions in particular: a) a non-transiting internal perturber whose interaction 
with the transiting planet is negligible, b) a non-transiting external perturber that is both 
eccentric and very well-separated (with a semimajor axis ratio > 5) from the transiting 
planet, c) the secular case with both planets on nearly circular orbits, and d) the resonant 
case with both planets on circular orbits. In all cases, the authors analytically estimate the 
timing variations caused by the non-transiting planet. For the middle two cases, the agree- 
ment between the theory and N-body simulations is particularly promising (their Figs. 2-3). 
The authors thereby help establish a collection of quantitative estimates for TTV signal 
amplitudes and provide details of their derivations. 



Nesvorny fc Morbidelli (2008), Nesvorny (2009) and Nesvorny fc Beauge (2010) devel- 



oped semianalytic expressions to approximate the magnitude of TTV signals from planets 
with given masses and orbital parameters. They also attempt to solve the inverse problem in 
a manner that is ~ 10^ faster to compute than direct N-body simulations in the planar case 
with a circular transiting planet and non-transiting external planet. [Nesvorny fc Morbidelli 



(2008) provide a formulation for solving the inverse problem in low and mid-eccentricity 
0.4) secular regimes (their Fig. 4) given a "cutoff parameter" which determine the num- 



ber of Fourier terms to use in their analysis. Nesvorny (2009) extended these results to the 



case of eccentric transiting planets and non-coplanar planets, and for these systems illus- 
trated the phase space regions in which he can solve the inverse problem to various levels 
of accuracy (his Figs. 1 and 3). Such expressions provide a valuable guide for determining 
which configurations should be explained in more detail although the models have trouble 
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reproducing the correct behavior at strong MMRs. 



Several studies have addressed other important aspects of TTVs. Ford & Holman (2007) 



considered how TTVs can be used to determine the existence of and characterize Trojan 
planets - planets that approximately share the same orbit as an observed transiting planet. 
They analytically estimated the resulting TTV signal and demonstrate the goodness of their 
relation with N-body simulations. A related system architecture, which, in principle, could 
be identified through TTVs, is that of a transiting planet containing a large (terrestrial-mass 
or larger) satellite. During their orbits, the Earth leads or trails the Earth-Moon barycenter 
by up to 2.5 minutes, and Saturn leads or trails the Saturn-Titan barycenter by up to 30 



seconds (Deeg 2002). Hence, Deeg (2002) claimed that a (now realistic; see Knutson et 



al. 2007) detectabihty threshold of 10 s would be needed to reliably detect a moon in the 



Saturn-Titan system through TTVs. Simon et al. (2007) furthered the theory of detecting 



exomoons with TTVs by showing how a planet-satellite system can be represented by a 
single theoretical body on the planet-satellite line at the "photocenter" . Therefore, the 
phase space explored in this work may be extended to planet-moon systems through such 
a relation. Kipping (2009a) and Kipping (2009b) have described how exomoons can be 



characterized through the use of transit duration variations (TDVs) as a way of helping to 
break degeneracies of the inverse problem from TTVs alone. [Meschiari &: Laughlin (2010) 
recently demonstrated how Doppler velocity measurements may be combined with TTVs in 
order to help remove the degeneracies inherent in identifying the architecture of the system. 
They considered systems similar in configuration to HAT-P-13 and HD 40307, as well as one 
with a transiting giant planet and a terrestrial-mass companion trapped in low-order mean 
motion resonance. 



Both the high precision (a few seconds; Knutson et al.|20^ ) with which transits can now 
be measured and the close proximity (sometimes within 0.02 AU of their parent star ^^jpj ) 
of Hot Jupiters to their parent stars necessitate the consideration of physical effects (e.g. 
general relativistic precession, sunspots, internal gravity waves) that have previously been 
neglected in exoplanet studies. For TTVs, most of the effects cause changes which are either 
well below the current detectabihty threshold or act on timescales much greater than a few 



years. For example, Miralda-Escude (2002) finds that over the course of a typical observing 



campaign (several years), both general relativity and the quadrupole moment of a star will 
induce a pericenter change of a Hot Jupiter of much less than one degree. However, |Heyl| 



& Gladman (2007) propose a method for using TTVs in order to determine the J2 moment 
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of a parent star, and Pal & Kocsis (2008) focus on detecting the effect of general relativity 
from both TTV and TDV observations. 



1.4. Our Methodology 

This paper primarily explores TTV phase space regions which are likely to be obser- 
vationally relevant and are not well described by analytical theories. This choice dictates 
that we rely on N-body simulations. Running ~ 10'' simulations is not prohibitively time- 
consuming since each simulation runs for just < 10 yr. This regime includes a circular 
transiting Hot Jupiter and an exterior coplanar lower-mass planet, with no additional restric- 
tions. We emphasize the sensitivity of TTV profiles to the initial orbital angle configuration, 
an important relation that has often been glossed over in the literature. We also showcase 
the difficulty in providing any constraints on systems with under 50 transit observations, 
despite several investigators' best efforts. Finally, we caution investigators on the dangers 
of confusing two planets which are locked "in resonance" with two planets whose periods 
are roughly commensurate. These results expand upon the preliminary investigations from 



Veras & Ford (2009) 



In Section 2, we a) describe the simulations which serve as a template from which other 
plots in the paper may be compared, b) state our definitions, and c) provide qualitative 
estimates for the magnitudes involved in solving the inverse TTV problem. In Sections 3-6, 
we demonstrate the dependence of TTV signal amplitude on number of transit observations, 
orbital angle configuration, planetary masses, and time evolution of orbital elements, respec- 
tively. We introduce a method for modeling the shape of a TTV curve in Section 7, discuss 
radial velocity follow-up, light travel time, and extensions to this study in Section 8, and 
conclude in Section 9. The Appendix details the derivation of libration widths for some of 
the relevant resonances considered here. 



2. Fiducial Simulations 

As a first step towards understanding the large parameter space for 2-planet systems, we 
establish a fiducial set of simulations which broadly characterize the TTV signals produced 
from a transiting hot Jupiter and an external terrestrial-mass planet. We consider a Mj = 
IMj inner planet which transits its parent = IMq star on an initially circular orbit 
(cj = 0) at ttj = 0.05 AU, and an external Mg = IM^ planet with semimajor axis, Oq, and 
eccentricity, Cq. The subscripts "i" and "o" are abbreviations for "inner" and "outer". The 
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mean longitude, mean anomaly and longitude of pericenter will be denoted by A, 11 and w, 
respectively. For any system, only one of A or 11 needs to be specified; either determines a 
planet's location along a Keplerian orbit, and the former is often used in planar resonance 

. Henceforth, all stated orbital parameters are assumed to 
be initial values unless an explicit time dependence is included. 

After establishing values for aoj ai^to^Wo^ Ilj or Aj, and Ho or Aq, we integrate the 
system over 10 years and tabulate the times of transits, and subtract the best fit Keplerian 
model to generate a TTV curve. The deviations of the transit times from strict periodicity 
tend to form periodic patterns which can have timescales extending to the duration of the 
simulations. We display a sample of these patterns and periodicities in Fig. [TJover the course 
of the nominal Kepler mission lifetime (3.5 yrs) in order to motivate the other figures in this 
paper. The dots represent transit times, and the amplitude of their variations ranges over 
three orders of magnitude across the panels. The variety of the patterns seen hint at the 
difficulties in deriving masses and orbital parameters of the external planet from TTV curves 
alone. The top panel curve (ao/ctj = 3.6593, Cq = 0.596, Ho = 70°, Ilj = 0°,tz7o = 180°) may 
be fitted well with a single sinusoid. The curve in the next panel below (oo/aj = 2.3313, Cq = 
0.395, rij = Ho = 0°,Ci7o = 0°) exhibits several modulations in addition to a simple sinusoid. 
Up until ~ 1 yr, and between 2.5—3.5 yr, this curve appears to exhibit a long period (> 3.5 
yr) trend with a modulated sinusoid of amplitude less than 10 s. In between, at ~ 1 — 2.5 yr, 
the curve appears to oscillate with an amplitude of 20 s without modulation. The curve in 
the second panel from bottom (ao/oj = 1.5812, Cq = 0.204, Ilj = Ho = 0°,-a7o = 0°) appears 
to exhibit a very long periodic trend (> 3.5 yr) with a high amplitude (hundreds of seconds) 
and a sawtooth-like modulation of amplitude < 100 s . Note the distribution of spacings 
between the transit times on the modulated sinusoids. These spacings may provide hints as 
to the architecture of the system (see Section 7). The bottom panel (ao/c^i = 1.5812, Cq = 
0.219, Ilj = Ho = 0°,tt7o = 0°) displays TTVs which vary by thousands of seconds over the 
first ~ 2 yrs before abruptly taking on the form of a simple sinusoid with amplitude < 10^ s. 
This panel especially highlights the importance of considering a sufficient number of observed 
transits with any type of analysis in order to avoid spurious conclusions. We emphasize that 
the four curves in Fig. [T] are not representative of all curves henceforth studied, but rather 
illustrate different types of trends and patterns that one might uncover. 

In principle, we may select any combination of transits to mimic observed data. For 
our fiducial case, we select the TTV curves based on the first 10 years of transits (A^ ~ 874, 
where iV is defined as the number of transits). Figure [2] displays the results of integrations of 
120,000 systems. We sampled 400 logarithmically-spaced values of the semimajor axis ratio 
and 60 uniformly-spaced values of Co for each semimajor axis ratio. For every given pair 
of aoj ai and Cq values, we simulated 5 systems initialized with random mean anomalies of 



studies (Murray fc Dermott|[2000 
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Fig. 1. — A sample of transit timing variation (TTV) curves for four different systems over 
the course of the nominal Kepler mission lifetime (3.5 yr). Each dot represents a transit, and 
the vertical axis labels are shghtly offset for clarity, but do reflect the range of TTVs. The 
curves from top to bottom correspond to {ao/ai = 3.6593, — 0.596, Ho — 70°, tUo = 180°}, 
{tto/a, = 2.3313, Co = 0.395, = 0°, tUo = 0°}, {ao/ai = 1.5812, = 0.204, = 0°, = 0°}, 
and {ao/ai = 1.5812, = 0.219, n„ = O°,zuo = 0°}, with = 0° for all curves. Note the 
wide variety of patterns exhibited by the TTV curves. 



both planets and a random longitude of pericenter of the outer planet. The reported values 
represent the median root-mean-square (RMS) TTV deviation amplitude (in seconds) of each 
of these sets of 5 simulations. We only sampled systems guaranteed to be stable according 
to the Hill Stability Limit ( Gladmaii||1993 ) , i.e. where Co < en- The value of en is a function 
of both planets' masses, semimajor axes and eccentricities, and for Co < en, the orbits are 
guaranteed to never cross. As we vary these parameters, the bounding curve (where Co = en) 
on the figure will change slightly. We caution that i) some (unsampled) systems which lie 
above the curve in Fig. 1 may be stable, and that ii) some systems close to the Hill stability 
boundary which we do sample might be Lagrange unstable (where the outer planet drifts 
outward, generally causing an increase in the RMS TTV signal). We explore the extent of 
the systems that are likely to be Lagrange unstable in Section 8.1, but note briefly here that 
such systems are very unlikely to be observed. The inner semimajor axis ratio bound of 1.3 
was chosen to roughly correspond to the point where the Hill Stability curve intersects the 
X-axis. The outer semimajor axis ratio bound of 5.0 is arbitrarily chosen to allow one to 
consider highly hierarchical (widely-separated) systems. Additionally, 5.0 is the ratio beyond 
which analytic formulae achieve success at reproducing the RMS TTV amplitude given the 



orbital parameters of the exterior planet (Fig. 2 of Agol et al. 2005). 



Figure [2] is contoured on a logarithmic scale. If one takes 10s as the current detectability 
threshold for TTV signals, then any region on the plot that is not white, pink or red should 
contain a detectable signal. Hence, widely separated {ao/ai > 3) planets containing an 
external planet on a low-eccentricity (ci < 0.2) orbit produce currently undetectable signals. 
Additionally, for almost any Oo/ < 5 with an external planet whose orbital eccentricity is 
close to the Hill Stability limit, the resulting TTV signal will be easily (> 10^ s) detectable. 



The figure demonstrates that the signal amplitude is highly (on the 10 ^ 
sensitive to the semimajor axis ratio. 



10-2 AU scale) 



The fiame-like features on the plot, similar to those found in Fig. 5 of Agol et al. (2005) 



indicate regions where the two planets are near a period commensurability (PC) that can 
be expressed as a ratio between two small integers {p:q, where p,q < 20). The longest of 
these "flames of resonance" correspond to period commensurabilities of the form p:l, where 
p < 11. These PC locations approximate well regions where both planets might be locked in 
or reside just outside of a mean motion resonance (MMR). Whether or not the planets with 
a given p : q ratio are in a MMR is subject to individual system study, and is often highly 
dependent on its orbital angle architecture. In the exoplanet dynamics community, a widely 
used criterion for determining whether a system is "in" MMR is to determine if at least one 
"resonant angle" is librating. Consideration of multiple resonant angles for a given p:q can 



be found in Fig. 3 of Laughlin & Chambers (2001), Section 4.2 of Kley et al. (2005), Section 



4.1.1 of Raymond et al. (2008), Fig. 9 of Crida et al. (2008), and Section 5 of Fabrycky 
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Fig. 2. — "The Flames of Resonance"; the median RMS TTV amphtude (5'(Q)) for 5 
different initial orbital configurations and ~ 10 yrs-worth (A^ = 874) of consecutive transits 
for a transiting IMj hot Jupiter on a circular orbit at 0.05 AU and a IM^ external perturber 
with the orbital parameters indicated on the plot. The contour levels in seconds are: pink 
(0-1), red (1-3), orange (3-10), yellow (10-30), light green (30-100), olive (100-300), 
blue (300 — 1000), and purple (> 1000). Overplotted are resonant libration widths for selected 
PCs (Period Commensurabilities) . 



fc Murray- clay] ( |2010D . However, |Michtchenko et al.| ( |2008a| |b| claim that for a given 

a MMR is characterized just by one resonant angle, whereas in the planar case, the other 



independent angular variable should be a secular angle. Mardling (2008, 2010) asserts that 



for a given p:q, the resonant angle can be expressed in terms of a single index that traces 
the contribution from a particular semimajor axis ratio order in the gravitational potential. 

For our purposes, we can neglect the contribution of the inner planet's longitude of 
pericenter (wi) because = and the mass ratio of the outer planet to inner planet {Mo/Mi) 
is negligible (~ 0.001). Hence, ei{t) will remain low (< 0.01) as t increases, and the only 
resonant angle we consider is: 



^p,q 



p\o - qXi -{p- <l)^o- (1) 



The libration width of this resonance, defined here to be a range of Oq at a given Cq, 
bounds a region of phase space where this angle might librate. We use the formalism in 



Murray & Dermott (2000) to compute libration widths for a sampling of lst-4th order reso- 
nances, where order = p — q. Stan Dermott (2010, private communication) also provided us 
with the coefficients needed to compute libration widths for 5th-8th order resonant libration 
widths, with q = 1. The horizontal spacing between the thin vertical lines extending from 
the X-axis on the plots in Fig. [2] represent these libration widths. Notice that the spacing is 
negligible until Co > 6^/2. Further, the lines do not extend to the Hill Stabihty boundary. 



because of the Sundman convergence criterion ( Ferraz-Mello 1994 Sidlichovsky & Nesvorny 



1994). See the Appendix for additional details and a derivation of this libration width. 



The location and extent of the libration widths help confirm that pronounced RMS 
TTV amplitudes (= ^(Q) = S{t,ao/ai,eo,Ilo,^i,^o)) are likely caused by PCs and that 
MMR configurations are achieved only for a high enough eccentricity of the outer planet. 
Additionally, the libration widths for the 9:1 MMR demonstrate that two planets may be 
locked in that resonance despite the weakness of an 8th-order resonance (resonant strength 
scales as ~ e^""^). Such planets would likely produce a detectable TTV signal. Note also 
that the length of each "fiame" corresponds to the value of q. Planets locked in a high-order 
(> 10) MMR with g = 1 or 2 can produce a detectable signal (5'(Q) > 10s) for sufficiently 
large Cq. Although the highest-order MMR which would produce a detectable signal is highly 
dependent on the characteristics of the systems studied, planets in MMR of order 20 — 30 can 
in principle produce detectable signals if Cq is high enough. Thus, a transiting Hot Jupiter 
at 0.05 AU could have detectable TTVs due to an Earth-mass planet in the habitable zone 
of a K or M-type star. Because S'(Q) is highly (on the 10^'^ — 10^^ AU scale) sensitive to 
the semimajor axis ratio, detectable MMRs may help identify planets with particular orbital 
parameters. However, when considering high q MMRs, one should also take into account 



the time dependence of resonant locations; as a system evolves, libration widths would likely 
oscillate in ao(t) — Coit) space. 



Figure |2] displays the median root-mean-square (RMS) TTV deviation amplitude (in 
seconds) of five different orbital angle configurations whereas Fig. [s] reports the minimum 
and maximum deviations. We emphasize that the following conclusions are for our fiducial 
case, with a given perturbing mass of IMq. The differences in the two plots in Fig. |3] 
are striking. For a given Oq and Co, 5'(Q) may vary by orders of magnitude based on the 
initial values of the mean longitudes/mean anomalies and the outer planet longitude of 
pericenter [wi is meaningless when Cj = 0). The range of TTV signal amplitudes achieved 
from the different combinations of initial angles is comparable to the lower panel of Fig. [3| 
the maximum signal achieved. Further, the median signals are on average greater than the 
mean signals (not shown). The secular regime, away from PC, demonstrates little variation 
(at most a few seconds) in signal with initial orbital configuration. However, detections in 
these regimes are unlikely without the advent of better-precision observational techniques. 
The upper panel in Fig. |3| which displays the minimum signal achieved, best shows that for 
many strong (order < 7) MMRs with g = 1 or 2, the RMS signal is lower than that from 
the high-eccentricity (at en) secular regime. Therefore, for observations of a Hot Jupiter, 

1) moderate (tens of seconds) values of ^(Q) could indicate that Co may be as high as en, 

2) high amplitude (> 10'^) signals are not necessarily indicative of close proximity (within a 
few libration widths) to a strong PC, and 3) high amplitude (> 10^) S{Q,) values are always 
indicative of a high value (> O.Sen) of Cq. 

The choice of contour levels is arbitrary, but meant to show phase space structure and 
establish the detectability threshold of 10s. Increasing each contour level by one order of 
magnitude indicates that little structure can be discerned in the highest signal (purple) 
regions of Figs. [2] and [s] and that S{Q,) can reach well over lO^s. Contrastingly, lowering 
each contour level by one order of magnitude demonstrates clear structure in the lowest 
signal regime (pink) of Figs. |2] and [3] TTV curves in this regime have been successfully 



correlated with the mass and orbital elements of a hypothetical planet (Agol et al. 2005 



Nesvorny & Morbidelli 2008). However, characteristic exoplanet detection thresholds (~ 10 



s) exceed typical signal amplitudes in this regime. 



3. Correlations with the Number of Observations 



TTV signal amplitudes, 5'(Q), crucially depend on the number of observed transits, 
and the sampling rate. Although one can define TTVs solely as the difference between 
consecutive transit times, Agol et al. (2005) define TTVs as the deviation from an overall 
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RMS TTV (seconds) 




Semimajor Axis Ratio 




Fig. 3. — The minimum (upper panel) and maximum (lower panel) RMS TTV amplitude 
(^(Q)) from the 5 different initial orbital configurations sampled in Fig. |2} Note that (^(Q)) 
can vary by orders of magnitude depending on whether the minimum or maximum value is 
sampled. 
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linear fit to the transit times of a system. Tliis definition does not rely on having succes- 
sive transits. Before missions such as CoRoT and Kepler, ground-based observations have 
struggled to yield > 100. These already-operating space-based missions, however, hope to 
observe consecutive transits for a half year (CoRoT) or multiple years (Kepler). Therefore, 
we consider our fiducial system (Fig. |2| we henceforth use the median value when averaging 
over multiple random angle initial configurations) for various durations corresponding to 
= 874, 313, 100, 50, 30, and 10 (Figs. |4||6| and by considering consecutive transits. The 
first two values of listed correspond to observing campaigns of approximately 10 yrs and 
3.5 yrs (roughly the maximum and nominal lifetime of Kepler). 

Figure |4] illustrates the TTV amplitude variation as a function of A^. For the A^ = 10 
plot, almost no systems are detectable. At A^ = 30, only a few systems featuring high eccen- 
tricity may be detectable. For A^ > 100, systems close to a g = 1 PC generate signals which 
change little with A^ for Co < O.TbcH- However, S{Q) continues to increase for increasing A^ 
for Co values closer to ch- 

Because S'(Q) is sensitive to the semimajor axis ratio on scales potentially smaller 
than 10^'^ AU, higher resolution sampling of phase space will reveal additional features. 
Similarly, as A^ increases, larger sections of ao-Co phase space will exhibit TTV curves with 
clearly discernible periods. Hence, we have performed high resolution (200 points in x and 
120 points in y) explorations of the phase space around the 3:1 PC in Fig. [s] The figure 
demonstrates that as A^ increases, the contour outlines become sharper and other weaker 
commensurabilities which neighbor the 3:1 PC become more apparent. Additionally, the 
double- lobed feature on the 3:1 PC location at Co = 0.1 disappears with increasing A^. 

Importantly, the figure demonstrates that planets in a PC, and possibly in a MMR, do 
produce a distinct TTV signature, but not necessarily a high amplitude signature compared 
to its near-PC surroundings. At A^ = 30, ^(Q) at PC (yellow) is lower than the signal 
near PC (light green). At A^ = 313, S'(Q) at PC (light green) is lower than the signal near 
PC (olive), and this result is insensitive to further increases of A^. The near-PC regime 
is sharply divided from the secular regime (note the thin light green contour between the 
regimes). Additionally, after 10 yr, there is a sharp boundary (marked by sparse blue dots) 
between the near-PC regime and the high-eccentricity (at ~ ch) secular regime. Ten years 
might not be long enough for the some TTVs to produce periodic signals in the blue or 
purple areas. 

The variation of TTV signal with A^ is not necessarily monotonic, and is a function of 
the masses and orbital parameters of the planets. We can look more closely at this S{Q) vs. 
A^ dependence by fixing all other orbital parameters in specific cases. Although identifying 
"representative" systems is difficult because of the large size of the phase space, we sample 
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Fig. 4. — The Flames of Resonance (Fig. |2| for different numbers of consecutive transit 
observations; N = 874 (corresponding to ~ 10 yr), = 313 (corresponding to ~ 3.5 yr; 
the nominal Kepler lifetime), and = 100,50,30 and 10. The contour levels are the same 
as in Figs. [2] and [3] Overplotted are resonant libration widths for selected PCs (Period 
Commensurabilities). Note that for detectable signals (> 10 s), often several years-worth of 
observations are needed to pinpoint the true TTV amplitude. 
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Fig. 5. — High-resolution (200 points in x and 120 points in y) RMS TTV signal amplitudes 
for the region around the 3:1 PC for N = 874 (upper left), = 313 (upper middle), = 100 
(upper right), = 50 (lower left), = 30 (lower middle), and A^ = 10 (lower right). The 
contour levels are the same as in Figs. [2]and[3j Note that after several years of observations, 
near-PC signal amplitudes are higher than in-PC amplitudes. 
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8 different particular systems in each of three regimes: "in-PC", "near-PC" and "secular". 
Table [l] contains the ao/oj and Cq values for all 24 systems, and Fig. [6] plots S'(Q) vs for 
these systems. 

We first consider the "in-PC" systems plotted in the upper plot of Fig. [6j In order to 
model the signal amplitude of any in-PC system with continuous observational coverage, one 
must observe at least 50 transits. Only for > 50 do any of the curves level off. Generally, 
the higher the amplitude of >S'(Q), the greater the number of observations needed to level 
the curves (for a given perturber mass). Only the red (4:1), salmon (11:1) and purple (11:3) 
curves are level at iV = 874, and only for the blue (2:1) curve is S'(Q) decreasing. The 
signal curves are generally non-monotonic with N; additional observations may cause the 
signal amplitude to increase or decrease. Note additionally that although eccentricities are 
high {0A7eH < Co < O.SOen), the signal amplitudes are relatively low (5'(Q) < 200s) for 
all resonant curves for all t < 10 yr. Systems in a PC do not necessarily achieve signal 
amplitudes as high as those in other regimes, and preferentially have lower amplitudes, even 
after 10 yrs of observational sampling. 

The prospects for detecting systems near but not on PC is more promising. The middle 
panel of Fig. [6] indicates that most of the curves level off after !^ 300 (a little over 3 yrs). 
However, at = 800, one can see the aqua (5:2) and olive (3:1) curves beginning to trend 
upwards. The system near the 2:1 PC (blue) exhibits similar oscillatory behavior from the 
2:1 in-PC system in the upper panel. The two curves corresponding to the highest Cq values 
(salmon with Cq = O.Sden and magenta with Co = 0.806^) maintain S'(Q) < 20s for all A^. 
The salmon curves in the upper two panels of Fig. |6] contain the same value of Co and both 
level at 20s. However, the in-PC system levels off A^ ~ 230 whereas the near-PC system 
levels off at A^ ~ 110, reflecting the general trend of near-PC systems having robust values 
of S'(Q) with fewer observations. 

The two secular systems with the greatest signals in the lower panel of Fig. [6] (olive and 
magenta) undergo variations on the order of their initial values for at least 5 years. Both of 
these systems feature Cq > 0.80ei^ and Oo/aj < 2.2. The strongly hierarchical (flo/ctj = 4.405) 
system indicated by the gray curve raises its TTV signal by nearly one order of magnitude 
between A^ = 250 — 300, and then doubles this value over the next 6 years. The variation 
with A^ of all three of these systems dwarfs the variation by the other 5 systems on the 
plot, including 2 with high eccentricities (aqua and purple, with Cq = 0.73eH and O.Tlen 
respectively). These two high eccentricity systems show little (a few percent) variation for 
A^ > 200. Additionally, the signal for two of the lower eccentricity systems, indicated by 
the blue and red curves (co = OAAen and 0.25e//) appear to show negligible variation with 
A^ for 50 < A^ < 400. However, both signals begin to increase for A^ > 400 such that 
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Table 1. Sample Systems in 3 Regimes 
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Note. — 24 systems selected for individual analysis throughout the 
paper, with designations in Column 1. Column 2 lists the regime the 
system resides in (where PC is Period Commensurability) , and Col- 
umn 3 lists the corresponding PC, if applicable. The other columns 
list the initial semimajor axis ratio of the planets, the outer planet's 
eccentricity, and the color of the corresponding curve in Fig. |6j 
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N 

Fig. 6. — The RMS TTV signal amplitude as a function of number of consecutive transit 
observations (A^) for 8 systems in each of the following regimes: upper plot in-PC (Period 
Commensurability) , middle plot: near-PC, and lower plot secular regimes. The initial 
orbital angles of these systems are fixed; their other initial conditions are provided in Table 
[l] Note that generally TTV signal amplitude is a nonmonotonic function of A^. 
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S{N = 874) /S{N = 400) > 1.5. The qualitative difference in all these secular curves help 
illustrate that ^(Q) is non-trivially dependent on N, even for secular systems far from a PC. 

4. Correlations with Orbital Angles 

Although we have already shown in Fig. |2] and Section 2 how randomly chosen sets 
of initial orbital angles (mean anomalies and longitudes of pericenter) for a given and Co 
can produce order of magnitude variations in TTV signals, here we study this dependency 
systematically. In the following analyses, we keep all orbital parameters but one fixed in 
order to gain insight into the physical origin of the TTV signals. 

Although such insight is difficult to discern when considering how the entire (oo, Co) 
phase space varies with individual orbital angles, line plots for fixed (oq, Cq) values show 
more structure. First, however, we consider "flames" plots as a function of Uq- In each plot 
of Fig. [7[ the initial orbital angles are Ilj = = 0°, with II^ = 0°, 10°, ...80°. Systems with 
planets initially at conjunction have high signal amplitude at high (near en) eccentricity, 
and, unlike in all other plots, saturate PCs at these eccentricities with S'(Q) > 10^ s at 
N = 874. For IIq = 10° — 80°, in each plot, several PCs produce signals which are over 
two orders of magnitude lower than 10'^ s at high e,,. The "flames" which produce these 
effects appear at different sets of PCs in each plot, and have no immediately recognizable 
pattern. The lack of an apparent pattern in these plots, despite fixing all but one parameter 
and despite the high value of A = 874, expresses well the difficulty in trying to use transit 
timing variations to deduce the properties of unseen planets. 

An alternative and perhaps more demonstrative method of exhibiting this signal de- 
crease at high eccentricities is through sequences of individual transit curves. Fig. [8] displays 
9 sequences (panels) of scaled curves, where each curve traces 3.5 yr, as Co is increased from 
to en at the 7:1 PC {ao/ai = 3.66). In each panel, Cq is increased from the top left curve 
to the bottom left, to the upper right to the bottom right. The numbers accompanying each 
curve represents the RMS TTV amplitude of that curve in seconds. The curves are scaled 
such that their extrema are the plot boundaries, so that one could better see detail and 
modulation. 

The figure demonstrates that at this PC: 1) although low eccentricity (co ~ 0) curves 
are highly structured and modulated, their signal amplitude is negligible, 2) at mid-to-high 
(co ~ O.Tch) eccentricity, the frequency of the curves suddenly increases, 3) in only some 
cases, eccentricities closest to en will exhibit long-period (> 3.5) yr signals of > 10^ s, 4) 
the transitions to these high signals is often sudden (caused by a difference of e < 0.018), 
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Fig. 7. — RMS TTV amplitude for Hi = Wo = 0° as Ho is increased in 10° increments from 
0° to 80° in the upper left, upper middle, upper right, middle left, center, middle right, lower 
left, lower middle and lower right panels, respectively. The contour levels are the same as in 
Figs. [2] and [3j All systems are sampled for = 874. Note the different locations in which 
"flames" arise in the blue and purple regions in the plots. 



- 22 - 



n, = o° 



Ho = w 



Uo = 20° 




= 60° 



n. = 70^ 






































































vz-v/v/V 












\/\7\/\7\/ 






























.00 3.50 
Time (yr) 



64 

58 

55 

52 

49 

47 

45 

45 

45 

25 

23 

23 

22 

23 

22 

23 

23 

23 

24 

24 

24 

25 

25 

26 

35 

2211 

5521 

2251 

5011 

136; 




= 80° 






















































































A ;oaV 








^A\^--OA 































4S 

50 

49 

46 

43 

42 

43 

47 

52 

58 

62 

63 

62 

61 

37 

34 

33 

35 

34 

34 

34 

34 

34 

35 

73 

135 

308 

358 

9731 

3421 



Fig. 8. — Transit curves over 3.5 yr close to the 7:1 PC for Ilj = = 0° as Ho is increased 
in 10° increments from 0° to 80° as indicated on top of the 9 panels. Within each panel, 
there are 60 systems sampled; Cq increases from in the top left of each panel downward, 
first in the left column and then the right column, until = en in the bottom right of each 
panel. The vertical range of each plot is scaled, and we provide the value of S{Q) beside 
each curve. Note the qualitatively different locations in eccentricity space of TTV signal 
forms in the right column of each panel. 
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and 5) in no case does the signal monotonically increase with increasing Cq. In particular, 
for Ho = 10°, 40°, and 60°, TTV curves exhibit low (< 10s) amplitudes at en- For the other 
values of Uq, TTV signals at these high eccentricities exhibit long (> 3.5 yr) period variations 
and very high (> lO'^s) amplitudes. This qualitative difference might vanish at high enough 
eccentricities, beyond the Hill Stability Limit. Note additionally that the eccentricity range 
which allows for high-frequency (featuring several crests over 3.5 yr) TTV signals varies 
significantly (over a factor of 2) for different IIo values. For IIo = 50°, the long-trend 
signals near en are punctured by an orderly, periodic 89 s signal amplitude, showcasing the 
unpredictability of TTV signals near MMR. 

The above analysis keeps Wo fixed. Although this angle refers to the outer planet's orbit 
and not the planet's location on that orbit, the value of uio may play a crucial role in the 
dynamics. In fact, systems with tUo = 90° and Wo = 180° (not plotted) do show completely 
different patterns of low amplitude (< 100s) "flames" near the Hill Stability Limit from 
those shown in Fig. [7j The starkest difference appears in the Ho = 0° case, where these 
flames are absent for Wo = 0°, fully saturate the Hill Stability region for Wo = 90°, and 
appear only for Oo/aj < 2.5 when zUo = 180°. 

The above analysis does not necessarily extend to other PCs. As the coefficients of 
the terms in disturbing function, or the different shapes of the libration widths in Fig. |2] 
would indicate, the resonant structure of each PC is different. Analysis of planets thought 
to be near or in a particular MMR hence would benefit from a high-resolution and complete 
exploration of the phase space at that location. 

In order to better determine if the variations with Ho in Fig. [7] have structure, we 
now consider the signal variations while fixing (ao,eo). We plot this variation for selected 
systems (I, II, VII, XVII, XX, XXII from Table Q in Fig. || where the sohd, dotted and 
dashed curves in each plot correspond to zUo = 0°, 90°, and 180°. The y-axis of all plots are 
logarithmic. At the coarse resolution of 10° per data point, the curves in the upper panels 
are vaguely oscillatory. For both the 2:1 (I) and 3:1 (II) systems, the S{Q) sharply (by 
several hundred seconds over 10°) peaks at Ho = 90°, but for a different value of Wo at each 
PC. At the 7:2 PC, this peak does not occur; curves for the three values of sampled have 
S'(Q) with = 874 values all within 10s of one another at Ho = 90°. All three upper panel 
plots indicate an anticorrelation for the blue and olive curves, with the most pronounced 
difference (over 400s) for the 3:1 PC at Ho = 90°. 

The variation of 5'(Q) with Ho and cCo for the secular systems sampled (lower panels 
of Fig. |9]) qualitatively differ from those for the in-PC systems. For the low eccentricity 
(co = 0.14e//) system XVII, the curves have fewer extrema, and have S'(Q) < 10s except 
for the dotted and dashed curves peaking at 30° and 40°. The other two secular systems 
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Fig. 9. — TTV RMS signal amplitude vs. Ho for systems at the 2:1, 3:1 and 7:2 PC (upper 
panels, left to right) and for secular systems XVII, XX and XXII from Table [l] (lower panels, 
left to right). In each plot, the solid, dotted and dashed curves correspond to zuo = 0°, 
zuo = 90° and zuo = 180°. All signal amplitudes were sampled after 10 yr of continuous 
observations (A^ = 874). Note the lack of a discernible correlation of Uq with RMS amplitude. 
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sampled are at high (cq > O.TSen) eccentricity, and both feature curves with multiple signal 
extrema, despite the 3-4 orders of magnitude difference in ^(Q) between those two plots. 
The greatest extrema occurs at 60° for the dotted curve in system XX and at 70° for the 
solid curve in system XXII. The apparent anticorrelation between the solid and dotted curves 
from the in-PC systems appears to vanish for these secular systems. 

In order to observe the variation of 5'(Q) with IIq at a higher resolution of Ho, we 
performed additional simulations. For the secular systems XVII, XVIII and XIX, we varied 
Ho over the entire [0°, 360°] range by sampling the angle at intervals of 0.1°, and report the 
results in Fig. [TOj The three panels show qualitatively different behavior. The bottom panel 
(with ^(Q) < 1 s) feature smooth curves, the top panel curves (with 2.8 s < S'(Q) < 4.4 s) 
show spikes and dips which could be hidden in broader resolution studies, and the middle 
panel (with S'(Q) up to 7 x 10'^ s) shows little structure whatsoever. The figure demonstrates 
that TTV amplitude can vary by seconds due to changes of a few degrees in Ho in some 
regimes, and vary by thousands of seconds due to changes of a fraction of a degree in Ho in 
others. 

We plot these curves primarily to demonstrate the difficulties inherent in solving for the 
mass and orbital parameters of the unseen planet. Because of the broad parameter space, one 
is hard-pressed to determine "representative" systems in each regime. Additionally, because 
of the semimajor axis resolution of our simulations, the typical difference between nominal 
MMR tto values and those from our simulations is ~ 10~^ — 10~^. This difference may be 
significant due to the sensitivity of TTV signal profiles on the planetary semimajor axis 
ratio. Thus, we conclude that any TTV signal will need individual analysis, as opposed to 
being characterized by a few summary statistics. Further, limits on the mass of the planets 
possible for a given TTV dataset must be mindful of the possibility that a putative planet 
could have orbital elements that result in a signal much smaller than is typical for a planet 
of a given mass, semimajor axis and eccentricity. 



5. Correlations with Masses 

Thus far, we have fixed the masses of the star and both planets in all systems studied. 
Because these three masses are hierarchical (Mq ^ Mj ^ Mq), varying any one of them by a 
factor of a few won't significantly alter the contour phase space structure of Fig. |2} However, 
mass variation might qualitatively affect individual systems at the edge of a secular, near-PC, 
or in-PC regime, or with an outer planet at a moderate-to-high (> O.ben) eccentricity. 

We consider four different external perturber masses (IM^, 5M^, lOM^, 5OM0) and five 
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Fig. 10. — High resolution plots of TTV RMS signal vs. Uq- Every 0.1° of Uq is sampled for 
secular systems XVII (upper panel), XVIII (middle panel) and XIX (lower panel). In each 
plot, the black, blue and green curves correspond to tUo = 0°, tOo = 90° and Wo = 180°. All 
signal amplitudes were sampled after 10 yr of continuous observations (A^ = 874). Note the 
wide variation in sensitivity of signal amplitude to Ho- 
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different transiting planet masses (O.lMj, 0.5Mj, IMj, 5Mj, lOMj). Because Hill Stability 
is a function of these masses, transiting inner planets more massive than IMj will most 
significantly restrict the semimajor axis-eccentricity space in which the system is guaranteed 
to be stable. Therefore, in Fig. 11, the phase space plotted is smaller than that of Fig. |2] 



for all 20 combinations of planetary masses. 

Increasing the mass of the outer planet increases the TTV signal in almost all areas 
of phase space. More subtly, increasing mass causes the signal amplitude at resonance to 
be higher than that from near-resonance, in contrast to Fig. [5] Fewer PCs signatures are 
discernible at higher-mass transiting planets; for Mi = lOMj, only PCs with g = 1 or 2 can 
be identified by inspection in the figure, independent of the value of M^. Correspondingly, 
one can distinguish the greatest number of PCs for the lowest Mj values. Additionally, for 
the most massive interior and exterior planets, the highest-signal stable (purple) region is 
greatest in extent. 

Unlike the correlation between TTV signal and orbital angles, the correlation between 
TTV signal and planetary mass is robustly exponential in most stable regimes. Plot- 
ted in Fig. |12| are 6 systems (I, IV, IX, XIV, XVII, XVIII) from Table [l| with the 



solid, dotted, dashed, dot-dashed and triple dot-dashed curves corresponding to Mj = 
O.lMj, 0.5Mj, IMj, 5Mj, and lOMj. All six systems illustrate explicitly a roughly expo- 
nential increase in ^(Q) as a function of Mq. The upper panel systems, which correspond 
to systems in the 2:1 and 8:1 PCs, demonstrate that ^(Q) decreases with increasing Mj, 
unlike in the panels below, and that for 8:1, all curves monotonically increase, whereas the 
curves for 2:1 do not. The middle panel of the figure displays 2:1 and 5:2 near-PC systems 
with Co = OAAch and Cq = OAOen, respectively. The difference in S{Q) upon varying the 
planetary masses is striking: 6 orders of magnitude for the 2:1 near-PC system, and just 2 
orders of magnitude for the 5:2 near-PC system. The secular systems in the lowest panel 
show two anomalous features: 1) for system XVII, S'(Q) for the Mi = lOMj curve is 2 orders 
of magnitude higher than any of the other Mi curves, 2) system XXIII joins the in-PC 8:1 
system as the only ones which contain a curve (solid; Mj = O.lMj) where S'(Q) decreases 
as M„ is increased from IMm to 50Mffi. 



6. Correlations with Orbital Parameter Evolution 

Over the course of a multi-year observing campaign, the semimajor axis and eccentricity 
of the transiting planet might vary noticeably. N-body simulations show that an exterior 
terrestrial-planet could induce semimajor axis and eccentricity variations of ~ 0.003 AU and 
~ 0.005, respectively, of the massive transiting Jupiter. If this variation is observed over 
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Fig. 11. — The median RMS TTV amplitude for 5 different initial orbital configurations for 
(horizontally) 1, 5, 10 and 50 external perturbers and (vertically) 0.1, 0.5, 1, 5 and 10 
Mj transiting planets. The contour levels are the same as in Figs. [2] and [3j All systems are 
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Fig. 12. — TTV RMS signal vs. Mo for systems at the 2:1 and 8:1 PC (upper pan- 
els, left to right), systems near the 2:1 and 5:2 PC (middle panels, left to right), and 
secular XVII and XXIII systems (lower panels, left to right) from Table [l] In each 
plot, the sohd, dotted, dashed, dot-dashed and triple dot-dashed curves correspond to 
Mi = 0.1Mj,0.5Mj, lMj,5Mj, and lOMj. All signal amplitudes represent the median val- 
ues from sampling 5 randomly chosen sets of initial orbital angles after 10 yr of continuous 
observations. Note that although there is a power- law- like dependence of RMS amplitude 
on Mo (unlike for Mj), this dependence breaks down for some orbital architectures. 
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time, then in principle these profiles could supplement TTVs as a method for identifying or 
confining the parameters of unseen planets. 



Away from PC, although various formulations (Veras & Armitage 2007) of Laplace 



Lagrange secular theory (see Murray & Dermott 2000) may be used to provide approxima- 
tions of Ciit) profiles which don't exceed low-to-moderate eccentricities ~ 0.2 — OAen, these 
secular timescales (often 10'^~^yr) exceed the duration of optimistic observing campaigns 
(10 yr) by many orders of magnitude. One should instead consider evolution on orbital 
timescales. One could apply perturbation theory without averaging over mean longitude in 
order to model these systems. Alternatively, suites of N-body simulations are prudent since 
the CPU time required to evolve two-planet systems on the timescale of 10 yr is on the order 

of /iS. 

At PC, and especially in MMR, theory may yield an additional useful constraint. Be- 
cause of the mass hierarchy in our system, and the initially circular orbit of Mj, the resonant 
state of the systems studied here is likely to be dominated by a single resonant angle of the 
form in Eq. ([T]). In this situation, constants of motion exist which relate just two orbital 
parameters to each other through conservation of angular momentum and energy. For any 



eccentricity-type of resonance, one of these constants, Ci, may be expressed (Veras 2007) in 
terms of aj(t) and ej(t) only (recall that Cj = ei{t = 0) and = aj(t = 0)): 



P 



l-e 



(p-q) 



(p-q) 



(2) 
(3) 



Another constant, C2, relates ai{t) to ao{t): 



Co = 



qMA M^M^ + Mo + Mi 



- \/ai 



(4) 



= y^^Y - (5) 



We rewrite Eqs. (l2])-([5]) as 
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Vl - ei{ty = 1+ J 

(6) 

which illustrates how a transiting planet's eccentricity is predicted to change as a function 
of the variation in the external perturber's semimajor axis for a given MMR. However, this 
expression does require knowledge of the initial semimajor axes and eccentricities of both 
planets. 

Assuming that the outer planet is much more massive than the inner planet, the unseen 
outer planet will experience greater variations of its orbital elements over time than the 
transiting inner planet. Here we consider only the semimajor axis ratio range of 1.3-3.2. For 
Co = 0, the variation of eo{t) typically does not exceed 0.01, but may reach 0.1 for systems 
near the 2:1 PC. The eo(t) profiles in this regime include sinusoidal, lightly modulated, 
and heavily modulated curves with characteristic periods ranging from a few days to a few 
years. For Co = O.Bch, typically 0.01 < eo{t) < 0.1, with a greater proportion of sinusoidal- 
like curves. For Co = 0.75eH, typically eo{t) > 0.1, with variations up to fa 0.6. Many eo{t) 
profiles have httle apparent structure, and are periodic only on long (> lOyr) timescales. The 
outer planets with the smallest eccentricity variation, at eo(t) ^ 0.01, have corresponding 
TTV RMS amplitudes of just 10s. RMS signal amphtudes of ~ lO^s may correspond to 
Co variations of < 0.10 or > 0.40. For Co = en, eo{t) profiles vary chaotically and have 
typical amplitudes > lO'^s. For Qo/ai > 3.2, eo{t) profiles are generally flat except for those 
corresponding to Cq ~ Ci^, in which they vary chaotically. 

Transiting systems at or near PCs is perhaps of the greatest observational interest. 
Therefore, we trace the time evolution of the resonant angle in each system to determine 
if and over what timescale, and with what amplitude, this angle librates. We consider the 
strongest (lowest-order) PCs for g = 1,2 and 3, and fix 11^ = zUg = 0° and Ug = 60°. For 
each PC, we sample the resonant argument evolution for all simulations with the three values 
of Qo/ai simulated which are closest to the value of {p/qf^^- Here, we describe the results 
qualitatively for this narrow region of resonant phase space. 

All these systems are near a PC, and at eccentricities near en are all within the hbration 
width of the corresponding commensur ability. However, only a fraction of these systems 
actually exhibit a hbrating resonant angle, and hence are "in" MMR. The resonant angle 
may librate for different timespans; throughout an observational campaign or periodically 
during that campaign. The importance of a system actually in MMR versus just having a 
period ratio close to that commensurability is arguable, but a topic we now address. 



Y 



\/ao{t) 



\ -1 
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First, we consider a commensurability that is saturated by resonant systems. We plot 



the systems sampled for the 2:1 PC in Fig. 13 The upper, middle and lower panels represent 



systems at the three different initial semimajor axis ratios (1.5812, 1.5865 and 1.5919), and 
the left and right panels represent the ranges > Co > 0.5eH and O.den > > en, with Co 
decreasing downward. Every row of each panel represents a different system; the left column 
represents the TTV signal, with corresponding amplitude in seconds to the left; the middle 
column represents the "apsidal angle" (= Wo — ^i); the right column represents the resonant 
angle (2Ao — Aj — Wo) in red and another angle (2Ao — Aj — Wi) in green. The middle and right 
columns have 0° — 180° plot ranges, whereas the left column's range is scaled according to 
S'(Q). The curves are shown for a time evolution of about 1.1 yr (corresponding to ~ 100). 

The red curve librates for all 180 systems in the figure, demonstrating that all these 
planets are locked in MMR, at least for about 1 yr. At the lowest values of Co, this libration 
occurs about 0°, and is hence said to be "symmetric" and "aligned". For other values of 
Co, the libration is about 45° and is hence said to be "antisymmetric" . The green and blue 
curves librate for just a few systems at each value. The libration of these angles means 



that the system is in Apsidal Corotation Resonance (ACR; Beauge & Michtchenko 2003 



Ferraz-Mello et al. 2003 Lee|2004; Kley et al. 


2005 


|Beauge et al. 


2006 


Michtchenko et al. 


2006a 


Voyatzis & Hadjidemetriou 


2006 


Sandor et al. 


2007 


Michtchenko et al. 


2008a|b). 



The location of these systems occurs at ~ 0.2ejj, but does vary sensitively on the value of 
Oq. The 2:1 MMR is the only resonance studied which features a librating apsidal angle. 

The primary difference between the upper, middle and lower panels are the amplitudes of 
the TTV signals, which can vary by several factors despite the small (< 0.006 AU) difference 
between values. The highest signals (with > 100s) appear in two eccentricity groupings 
in the upper and middle panels, and just one eccentricity grouping in the lower panels. The 
first group of high signals in the upper and middle panels ends when the libration of the 
resonant angle becomes antisymmetric and the TTV signal qualitatively changes character. 
The gap between groups of high signals for the upper panels is five times smaller than that 
for the middle panels. For 0.5 < Cq < en, the libration curves, TTV signal profiles, and TTV 
amplitudes are relatively insensitive (varying by a few percent) to changes in Oq, unlike for 
< Co < O.ben- In this low-eccentricity regime, the highest amplitudes seen in the middle 
panels occurs at ACR; the same is not true for the lower or upper panels. In none of those 
panels does the TTV signal profile noticeably change at ACR. 

Unlike for the 2:1 commensurability, other p:g commensurabilities feature systems that 
are not in resonance, a situation corresponding to a circulation of the red angle. For all 
q = 1, and p = 2, ...11, only the 6:1 PC contains no MMRs. Systems also appear in MMRs 
for the 3:2, 5:2, 7:2, 9:2, 7:3, 8:3, 10:3, and 11:3 PCs, but none for the 5:3 PC. However, recall 
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Fig. 13. — Transit curves, apsidal angle evolution, and resonant angle (0) evolution over 
1.1 yr {N ^ 100) close to the 2:1 PC for = ro^ = ro^ = 0° and = 60°. All x-axes 
are in time. The values of ao/cti for the upper, middle and lower panels are 1.5812, 1.5865 
and 1.5919, respectively. In a given pair of horizontal panels, Cq is incremented uniformly 
from to en top to bottom, left to right. The RMS amplitude is given in seconds besides 
each transit plot (whose range is appropriately scaled) , and the range for the angle plots are 
[0, 180°]. The red curve plots Eq. (|T|) and the green the other 2:1 PC angle. The blue curve 
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importantly that these results hold for a narrow region of phase space where Hi = Wo = 0° 
and Ho = 60° and for just one set of masses. 

As Co increases, the transition between non-resonant and resonant systems produces a 
sudden change in TTV signal amplitude (by a factor of at least several) for all commensura- 
bilities where this transition can be seen. Therefore, TTV signals could clearly distinguish 
resonant from non-resonant systems, even if both are within the libration width of a par- 
ticular commensurability. Further, the transition between symmetric and antisymmetric 
libration corresponds to a sudden (by a factor of at least 2) change in TTV signal in several 
cases around the 2:1, 3:1, 4:1, 5:1, 7:1, 8:1 and 10:1 PCs. The signal jump caused by this shift 
in libration center is typically not as great as the jump caused by the transition to resonance 
from circulation. For 0.75e// < Cq < en, only the 2:1, 5:1, 7:1, 3:2 and 9:2 PC systems are 
in MMR. At the other commensurabilities, the resonant angle varies in a chaotic fashion, a 
qualitatively different behavior from the smooth circulation (featuring a continuous curve) 
found at non-resonant lower-eccentricity regimes. 



7. Correlation with TTV Signal Shape 

In previous sections, we have focused our investigations on TTV signal amplitudes. 
Their importance stems from the ability of observers to distinguish a TTV curve from mea- 
surement uncertainty. As we have demonstrated, the TTV signal amplitude is highly sen- 
sitive to the orbital parameters of both planets, and is of limited utility when attempting 
to identify an unseen planet. In order to help break the degeneracy, one can employ the 
shape of the TTV curve. As seen in Fig. [T| TTV curves may take on a variety of forms. By 
incorporating this form into a TTV analysis, one may be able to better pinpoint the mass 
and orbital parameters of the unseen planet. 

A quantifiable way of obtaining TTV shape data is to consider the autocorrelation 
function ^ at a given (time) lag L: 



Al = — — ( ( ) 



fc=l 



X 



k 



where Xk represents the deviation of the kth transit time from the constant period model (to 
be distinguished from S'(Q) at = A;, which represents the RMS amplitude from the first 
kth observations). The function Al is bounded by [—1, 1] and can be computed solely from 
the transit observation data alone. In principle, Al provides constraints, in addition to 
S'(Q), on the mass and orbital parameters of the unseen planet. 



- 35 - 



We find that the autocorrelation function may provide a useful summary of the TTV 



signal shape (see Figs. 14 15). Figure 14 plots the autocorrelation as contours for 9 different 
values of L (2, 3, 4, 5, 6, 20, 50, 100, 200) for = 313 3.5 yrs) and Y[o = ni = Wo = 0°. 
One notes immediately that although some "flame" -like features from Fig. [2] are apparent for 
low (L < 6) and high (L = 200) lags, the contours generally sculpt the phase space differently 
than does 5'(Q). This difference is most apparent in the low eccentricity (e,, < 0.2) regimes, 
which demonstrate rich structure in autocorrelation space. For L = 100, sharp "flames" 
puncture this regime for initial semimajor axis ratios of up to ~ 4. This low-eccentricity 
regime is transformed drastically as one increases L from L = 2 to L = 6. Note importantly 
that a purple flume (with Al close to unity) is almost precisely centered on the L:l PC for 
L = 2,3,4,5, and 6 (corresponding to ao/oj = 1.59,2.08,2.52,2.92,3.30 and 3.66). Such a 
correlation is a promising sign of the utility of Al in helping to diagnose plausible integrations 
for future study. 

The above contours fixed the initial orbital angles, which we have shown to sensitively 
affect the TTV signal amplitude. When we instead fix the lag, e.g. at L = 3, and vary Uq 
as in Fig. [7| then at the contour resolution of Fig. 14, the change in initial angles appear 



to have no effect on the value of Al- However, at finer resolutions, there is a variation. For 
some systems, this variation in Al strongly mimics the variation in ^(Q) as Uq is changed. 
Fig. 15 displays three secular systems: the system (XVII from Table [T]) in the upper panel 



with L = 3 shows remarkable agreement in the shapes and relative extrema of the solid, 
dotted and dashed curves (corresponding to Wo = 0°, Wo = 90° and tUo = 180°) to the 
signal from the bottom left panel of Fig. [9j In contrast, the middle panels, which feature 
secular system XVIII, does not share the same agreement with the corresponding signal (not 
shown). In between these two extremes, the bottom panels, featuring secular system XIX, 
illustrate different levels of correlation to the signal variation (not shown): strong for the 
dashed curve, moderate for the dotted curve, and none for the solid curve. Therefore, the 
RMS amplitude and autocorrelation value might be strongly or weakly correlated depending 
on the particular orbital configuration. The right panels of Fig. [15] demonstrate how the 
shape and magnitude of the Al curves can qualitatively change for two different (arbitrarily 
chosen) values of L. 

One may consider how Al varies with planetary mass as well as initial orbital angles. 
In many cases, we do find a correlation between the variation of Al and S'(Q) with Mg, even 
if the agreement fails to come close to that from some orbital angle regimes. The variations 
induced by changing Mq for a given L tend to be greater by an order of magnitude than 
those induced from orbital angles alone. In a sense, the variation caused by the angles is a 
modulation to those caused by the masses. 
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Fig. 14. — The autocorrelation function Al at lags corresponding to L = 
2, 3, 4, 5, 6, 20, 50, 100 and 200 from left to right, top to bottom, for Hi = Uo = Wi = Wo = 0° 
and N = 313 (~ 3.5 yrs). The contours are pink (—1.0 to —0.75), red (—0.75 to —0.50), 
orange (-0.5 to -0.25), yellow (-0.25 to 0), light green (0 to 0.25), olive (0.25 to 0.5), 
blue (0.5 to 0.75), and purple (0.75 to 1.0). Note that each plot represents an independent 
measure on the TTV signal, and the first few lags trace out the locations of PCs. 
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Fig. 15. — Autocorrelation ^3 (left panels) and ^100 (right panels) for the three secular 
systems XVII (upper panels), XVIII (middle panels) and XIX (lower panels) from Table [T| 
In each plot, the solid, dotted and dashed curves correspond to Wo = 0° , zuo = 90° and 
Ci7o = 180°. Here, = 874. Note the similarity in form of the ^3 for system XVII and 
its RMS amplitude (bottom left curve of Fig. [o]), and that when L is increased to 100, Al 
changes both form and magnitude for most secular systems. 
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These preliminary results suggest a promising avenue in which to pursue the inverse 
TTV problem, particularly if additional effects discussed in Section |8] are coupled with am- 
plitude and shape data. Future work will entail identifying the critical lags which produce 
the autocorrelation extrema, quantifying how low {L < 10) lags track L : 1 PCs, and incor- 
porating S'(Q) values and Al for iV — 1 values of L all as independent parameters to help 
remove the degeneracy of the inverse TTV problem. 



8. Discussion 



8.1. Lagrange Unstable Systems 

The highest signature TTV amplitudes (^(Q) > 1000s), like those generated from the 
curve in the fourth panel of Fig. [T| might be the result of erratic changes in the orbital pa- 
rameters of the planets. First, we confirmed that the highest eccentricity systems simulated 
for this work are provably Hill stable. In principle, such systems might be manifestly unsta- 
ble even though they satisfy Hill Stability. Although the planetary orbits will never cross 
in a Hill Stable system, the outer planet might get ejected because the outer planet is not 
bounded, or is Lagrange Unstable. If this is the case, then the timescale for this ejection is 
expected to be orders of magnitude less than the lifetime of the system, so observing planets 
in such configurations would be rare. 

We attempt to place an upper limit on the fraction of systems which exhibit this mani- 
festly unstable behavior by investigating configurations which produce sudden time variations 
of orbital elements such as the semimajor axis, periastron and apastron, mean longitude and 
longitude of pericenter. Such estimates may be used as a Lagrange stability filter in addi- 
tional statistical studies of TTVs. We find that comparing the maximum semimajor axis 
difference or ratio of both planets over the length of the simulation with their initial values 
provides a useful diagnostic. In particular, we define 



_ max (ao(t) - ai{t)) 
X= • (8) 

Therefore we compute x for every one of the 120,000 simulations studied in our fiducial con- 
figuration of Figure. [2| plus in selected configurations from Figure 11 representing extreme 



values of the transiting and exterior planet masses. We report the percent of all systems 



with X foss than given values from 1.1 to 10.0 in the left panel of Figure [T6| The figure 
demonstrates that the large majority of Hill Stable systems have planets which do not vary 
their semimajor axis difference by more than 10%; in these relatively quiescent systems, the 
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outer planet is unlikely to be ejected. All but a few systems simulated with a transiting 
planet of O.lMj exhibit this quiescent behavior. Alternatively, systems with a more massive 
transiting planet are the most likely to exhibit instability. Nearly 20% of all systems sam- 
pled with a lOMj transiting planet exhibit changes in their semimajor axis ratio approaching 
1000%. Such systems could become unstable (in the sense of the outer planet being ejected) 
on timescales much shorter than the system lifetime. 

In order to help determine what Hill-stable initial conditions most likely generate this 
suggested Lagrange instability, one can consider the initial eccentricity of the outer planet. 
We discover that highly eccentric outer planets close to the Hill Stability limit produce the 
sudden changes in semimajor axis characterized by x ^ 2.0. In order to quantify this finding. 



we plot in the right panel of Figure 16 the median (dots) and mean (squares) values of Cq as 
a percent of ch which cause x to be less than fixed values for the same mass configurations 
as in the left panel. Note that these eccentricities correspond to the highest TTV signals 
('S'(Q) ^ 1000s; purple and blue regions in "flames" plots) in the fiducial case. We caution 
that these results are statistical in nature. The evolution of a particular systems is highly 
dependent on the initial orbital element configuration and will likely require long term (3> 10 
yr) study in order to determine whether they will individually tend towards instability. We 
find this sensitive dependency by noting that for Co = en and for given (initial) sets of ao/ai, 
altering initial sets of orbital angles can cause x to vary from i=s 1.1 to ~ 10.0. 



8.2. Radial Velocity Follow-Up 

Ideally, the existence of unseen planetary companions suggested from TTVs can be 
confirmed with radial-velocity observations. However, the low mass and orbital separation 
of such unseen perturbers may preclude them from radial- velocity detection (as well as 
observational limitations due to stellar properties). In this section, we estimate the maxi- 
mum magnitude of the radial-velocity semiamplitude produced by an external pertuber in a 
guaranteed-to-be-stable multi-planet exosystem with a hot Jupiter. Because radial velocity 
semiamplitude = i^o oc 1/ a/1 — e^, a planet could most easily be detected by radial velocity 



surveys when its eccentricity is highest. At the Hill Stability limit (Gladman 1993), 



a 



3'^/Vo/ii AtoA*i (11/^j + 7/Xo) 



ifJ'o + f^z 



^4/3 



3 (/io + /ii 



(9) 
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Fig. 16. — Possible tracers of Lagrange Instability in Hill Stable systems. The left panel 
displays the fraction of all 120,000 systems from Fig. |2] (solid black curves) as well as 



from the configurations in Fig. 11 corresponding to (Mj = 10Mj,Mo = IM^ ; blue 



dotted), {Mi = 10Mj,Mo = SOM® ; green dashed), {Mi = 0.1Mj,Mo = IM® ; red 
dot dashed) and (Mj = 0.1Mj,Mo = 50M® ; magenta triple-dot-dashed) whose value of 
X = max (ao(t) — ai{t)) /{ao — ai) is less than the values on the X-axis. The right panel 
shows similarly colored curves but for the median (dots) and mean (squares) values of Co/ch, 
expressed as a percent, at which x exceeds the value on the X-axis. The figure demonstrates 
that the fraction of Lagrange Unstable systems is likely to be at the few percent level for 
our fiducial simulations, and correspond to the purple and blue regions of the "flames" plots 
(Fig. i- 
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where a = ai/ao, fit = Mi/M^, and fio = Mo/M^^. Further, the radial velocity semiamplitude 



of the outer planet in the coplanar edge-on case is (Lee & Peale 2003): 



1/3 



(10) 



° GiM^ + Mi + MoY 



Therefore, the maximum semiamplitude of the external perturber, assuming ^ 



Mi, Mo, is 



GAL ,,2 



-1/2 



:i2i 



This limiting equation, plotted in Fig. 17, assumes that 1) the system is seen edge-on and 
2) Co = en- The figure, which displays curves for given masses of 1Mq),5Mi^, 1OM®,2OM0, 
SOM® and IMsatum (in ascending order of the curves on the plot), demonstrates that given a 
hot Jupiter at 0.05 AU around a Solar-like star, Ko is a stronger function of than of a. A 
single Earth-Mass external perturber is not yet detectable by radial velocity measurements, 
but a Super-Earth may be detectable by this means, independent of planetary separation. 
Once detected by radial velocity measurements, those observations can be combined with 



TTVs in order to help characterize the architecture of the system (Meschiari & Laughlin 
20T0|. 



8.3. Light Travel Time 

The variation in Light Travel Time (LTT) during each transit may exceed the previ- 
ously cited detectability threshold of 10 s depending on the planetary masses and orbital 
configuration. Therefore, TTV models might need to incorporate leading-order LTT effects 
in order to correctly describe the motion. All our simulations incorporate this effect. Here 
we show where the effect might be important. In order to help determine when LTT effects 
may be neglected, we plot this contribution in Fig. [18] Note the domain of the plot is 
three times as wide as those from Figs. |2] and [3} but still sampled with 400 points along 



the X-axis. Figure 18 quantifies the contribution of LTT effects for = 874 (left panels) 



and A^ = 50 (right panel). The contours plot the quantity log (S'(Q)/T(Q)), where T(Q) 
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Fig. 17. — The maximum semiamplitude of the radial velocity signature of external per- 
turbers of masses equaling 1M^,5M^, 1OM0,2OM0 and SOM® and IMsatum, corresponding 
to the solid/blue, dotted/green, dashed/red, dot-dashed/magenta, triple dot-dashed/orange 
and long dashed/aqua curves, respectively. Note that, in the TTV regime studied in this 
work, a single Earth-Mass external perturber is not yet detectable by radial velocity mea- 
surements, but a Super-Earth may be detectable by this means. 
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is the TTV RMS amplitude produced by LTT effects alone. The figure demonstrates that 
the LTT's contribution for = for an idealized 10 yrs of coverage can be at the few 
percent to tens of percent level, but does not dominate the TTV signal. However, for a more 
massive external perturber, LTT can dominate the signal, especially for hierarchical (widely 
separated) systems with an outer planet eccentricity up to O.den- 



8.4. Extensions 



This study has attempted to thoroughly quantify the TTV properties in an observation- 
ally relevant region of phase space. However, additional regions of interest exist, and future 
studies could help characterize alternate architectures. Motivated by the recent discoveries 



of retrograde transiting planets, Payne et al. (2010) complements this work by considering 



the nonzero relative inclinations between two planets, and with a transiting planet that can 
be either internal or external to the additional perturber in the system. Subsequent studies 
aim to characterize the prospects for constraining the frequency of multiple stellar systems 
among planet-hosting stars with TTVs ( |Montalt"o 2010), to detect Trojan perturbers (e.g. 
Haghighipour & Capen 2010, in preparation), to break the degeneracy in the inverse problem 



by observing either each of two planets transiting in the same system (Ragozzine & Holman 



2010) or a transiting planet and a transiting moon ( Kipping||2010b ), and to combine Markov 
Chain Monte Carlo analyses of radial velocity data with TTVs for the recently scrutinized 
system HAT-P-13 (e.g. Payne & Ford 2010, submitted). Further, as already mentioned in 
Section 1, Transit Duration Variations (TDVs) may help constrain the degeneracies from 



TTVs. Kipping (2010a) has recently presented new analytical formula for TDVs, and Dvo 



rak et al. (2010) demonstrates a correlation between TDVs and the relative inclination of the 



planets in another system of recent interest, CoRoT-7. Beyond these studies, other possibil- 
ities abound. TTVs in systems with three or more planets will provide both an opportunity 
and a challenge for theorists to explain and interpret, and potentially habitable planets might 
best be identified through TTVs orbiting K or M-type stars. Additionally, non-continuous 
transit observations, which are particularly relevant for ground-based observatories, may 
have different dependencies on signal amplitude and shape than those displayed here, and 
will likely be important to consider on a case-by-case basis. 



9. Conclusion 



Our goal was to illustrate and quantify the challenges involved in characterizing the mass 
and orbital parameters of an external perturber in a system with an observed transiting hot 
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Fig. 18. — Median log of the ratio of RMS amplitude to the contribution from Light Travel 
Time (LTT) for 5 different initial orbital configurations for 874 (left panels) and 50 (right 
panel) consecutive transits, and for an outer planet mass of one (upper panels) or ten (lower 
panel) Earth masses. Note that LTT makes a contribution of at least a few percent in 
most areas of phase space explored here, and dominates the most hierarchical and lowest 
eccentricity secular regime for Super-Earth external perturbers. 
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Jupiter through the use of Transit Timing Variations. TTV signals may vary by orders of 
magnitude due to ~ 10~^ AU shifts in semimajor axes, ~ 0.005 shifts in eccentricity, or 
~ 1° shifts in orbital angles. However, TTVs are a sensitive probe of hidden planets, and 
can suggest the existence of external perturbers due to signatures produced at high-order 
period commensurabilities (PCs; where p < 20 and ^ 3). More specific conclusions 
from this work are: 1) moderate (tens of seconds) values of signal amplitudes could indicate 
that the outer planet eccentricity could be as high as the Hill eccentricity, 2) high amplitude 
(> 10^ s) signals don't necessarily imply close proximity to a PC, 3) high amplitude (> 10^ 
s) signals are always indicative of a high value (at least 4/5ths of the Hill eccentricity) of the 
outer planet eccentricity, 4) near-PC amplitudes are often higher than in-PC amplitudes, a 
result largely independent of the number of observed transits, A^, 5) signal amplitude is a 
non-monotonic function of N, 6) increasing the external mass generally increases the TTV 
signal, 7) increasing the internal (hot Jupiter) mass generally decreases the number of easily 
discernible PCs, and 8) the distinctive signal amplitudes for systems at PCs don't necessarily 
imply that those systems are actually captured in a mean motion resonance. We propose 
using the shape data of a TTV curve through the autocorrelation function as a method to 
help characterize the external planets, and believe that the method could provide a promising 
avenue of future study. 

We thank the referee for insightful observations and Eric Agol, Dan Fabrycky, Nader 
Haghighipour, Matt Holman, David Nesvorny, and Jason Steffen for valuable and extensive 
discussions. This material is based upon work supported by the National Science Foundation 
under Grant No. 0707203. 
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A. Appendix 

Here we derive the estimates for the hbration widths drawn in Fig. [2] for selected PCs of 
up to 8th order and discuss the goodness of the approximation. We take "hbration width" 
to mean roughly the Oq range for a given Cq in which MMR locking between two planets is 
likely. Because of the hierarchical nature of the three masses for the systems we consider, 
and the initially circular orbit of the more massive planet, we assume that the one argument 
given by Eq. ([T| of the disturbing function will dominate the gravitational potential. This 
approximation is vital to obtaining a tractable analytic expression. 

We can assess the goodness of the approximation of neglecting terms that include ei{t) 



by using Veras's (2007) term-based integrator. We tested systems locked in the 3:1 and 5:1 
MMRs by including different numbers of terms up to 4th-order. For all systems, we set 
= Mq, Mi = Mj, Mo = Me, = 0.05, and = ai{p/qfl^. For the 3:1 MMR, we 
set Co = 0.30, \i = 100°, Ao = O",^^ = 340° and = 180° in order to produce resonant 
behavior. For the 5:1 MMR, we set Co = 0.45, = 5°, Ao = 170°, Wi = 3° and zuo = 19° in 
order to produce resonant behavior. Importantly, for the 3:1 MMR, we set Cj = 0.00, but 
for the 5:1 MMR, we set Cj = 0.01. Initializing the hot Jupiter with a non-zero eccentricity 
might cause additional disturbing function terms to become significant, and is realistic given 
that an external terrestrial-mass perturber can typically force the hot Jupiter's eccentricity 
to values of ^ 0.005. 

We test these systems by including all arguments that, in their coefficients, include 
ei{t) (case I), and by neglecting all such arguments (case II). For both cases and for both 
MMRs, the resonant angle librates, and does so about 180°. In Case I, for the 3:1 and 5:1 
MMRs, the hbration amplitudes are ~ 55° and ~ 65°; Case II changes these values by at 
most a few degrees. The range of eo{t) over 3 years in Cases I and II is 0.3040 — 0.2639 = 
0.0401 and 0.3035 - 0.2586 = 0.0449 for the 3:1 MMR, and 0.4669 - 0.4346 = 0.0323 and 
0.4669 — 0.4346 = 0.0323 for the 5:1 MMR. These results demonstrate that neglecting these 
additional arguments affects the motion in too small a manner to necessitate inclusion in 
this general, qualitative study. 

With a single disturbing function argument, 0, we can approximate a hbration width 



according to the same prescription found in Murray & Dermott (2000). The angle behaves 
like a pendulum so that: 

= w^sin0 (Al) 



One can find the value of co through a combination of 1) the definition of 0, 2) use of the 
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disturbing functions, R = R{(f)), and 3) use of Lagrange's Planetary Equations (Murray & 



Dermott 2000). We can define two disturbing functions for the systems we study here as: 



i?i = — er'' (a/d + /.)cos, 



R2 



a 
GMq 

an 



el-'ifd + fe) COS. 



(A2) 
(A3) 



where G is the Gravitational Constant, a = ai/ao, and fd, fid and fe are functions of a and 
Co. By defining ak = G {M^ + Mk), for k = i,o, we find: 



puo - qrii 
3 1/2^-5/2 



-3p dRi 



do 



' 1/2 -5/2. 

-ga/ a- a. 



"iq dRi 



dXo af dXi 



(A4) 
(A5) 

(A6) 



so that: 



[afd + fid] + [fd + fe 



(AT) 



The energy of a pendulum is = (1/2)0^ + 2ijj'^ sin^((/)/2). Hence, the maximum energy 
is 2a;^, and equating this value with E yields a relation between and (p. One can combine 
this relation with the Lagrange Planetary Equation for a in order to obtain: 



pGaJ^Mi 

- 1/2 ' 
ooaiao 



,p-g 



— fd + fii 



smi 



cos (0/2) ^ 

which can be integrated to finally obtain the (maximum) libration width 5a'- 



(A8) 



V3 ay\M,+Moy^^ 



Co 



(p-'/)/2 



Oifd + fid) 



afd + frd) + '-^{fd + fe) 



(A9) 



in the approximation Mo = 0, 
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(AlO) 



For most MMRs, fid = /e = 0. However, when q = I, one of the most important classes 
of MMRs for TTVs, both fid and fe are nonzero. Formulas for fd, fid and fe are provided 
in Murray & Dermott (2000) up to 4th order in Cq. Stan Dermott, private communication, 



has provided us with terms up to 8th order. Explicit formulas for fd in terms of a may be 



found in Veras & Armitage (2007) 



The eccentricity at which Eq. (AlO ) holds is restricted by the Sundman criterion (Ferraz 



Mello 1994 Sidlichovsky & Nesvorny 1994), a fundamental convergence criterion on the 



planar expansion of the disturbing function from Ellis & Murray (2000). This criterion can 
be expressed as: 



aiD{ei) < aod{eo), 



{All] 



where 



D{y) = \/l + y"^ cosh z + y + sinh z, 
d{y) = a/1 + cosh z — y — sinh z 



(A12) 
(A13) 



such that z is implicitly defined as z = g cosh 2;. This restriction prevents computation of 



the libration width for all eccentricities up to the Hill Stability limit. Figure 5 of Nesvorny 



& Morbidelli (2008) shows a detailed view of the difference between the Hill Stability limit 



and the Sundman convergence limit for several different mass ratios. 
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